Chapter10
Zhang Chen
8/8/2018
library(rethinking)
library(MASS)Easy
# 10E1
log(0.35/(1-0.35))## [1] -0.6190392# 10E2
exp(3.2)/(1+exp(3.2))## [1] 0.9608343# 10E3
OR <- exp(1.7)10E3
A log-odds of 1.7 corresponds to odds ratio of 5.4739474. That is 447.5% increase in odds.
10E4
When events are counted on different temporal or spatial scales.
Medium
10M1
In the aggregated form, the likelihood is Binomial, with n being the total number of trials. In the disaggregated form, the likelihood is Bernoulli, with n being 1.
10M2
When the predicotr increases by one unit, the order (i.e., log) of the expected value of count number in the Possion distribution increases by 1.7.
10M3
Because the linear combination of predictors may extend from negative infinity to positive infinity, while the probability of a certain event happens must be between 0 and 1. The logit link function allows the probability parameter to be modelled as a linear combination of predictors.
10M4
In Possion distribution we are modeling the expected value of count, which has to be positive. The log link function ensures that.
10M5
It would imply that the mean of a Possion GLM can only be between 0 and 1.
10M6
The binomial distribution has maximum entropy when only two events are possible, and the expected number of each event is assumed to be constant. The constrait for the possion distribution is essentially the same.
Hard
10H1
# trace plot for stan model
plot(m10.4_stan)
# compare precis output
precis(m10.4_map, depth=2)## Mean StdDev 5.5% 94.5%
## a[1] -0.73 0.27 -1.16 -0.30
## a[2] 6.67 3.61 0.90 12.45
## a[3] -1.03 0.28 -1.48 -0.59
## a[4] -1.03 0.28 -1.48 -0.59
## a[5] -0.73 0.27 -1.16 -0.30
## a[6] 0.21 0.27 -0.21 0.64
## a[7] 1.75 0.38 1.14 2.37
## bp 0.82 0.26 0.40 1.24
## bpC -0.13 0.30 -0.61 0.34precis(m10.4_stan, depth=2)## Warning in precis(m10.4_stan, depth = 2): There were 3 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a[1] -0.73 0.27 -1.15 -0.27 3692 1
## a[2] 11.03 5.28 3.42 18.34 2075 1
## a[3] -1.05 0.28 -1.49 -0.61 3091 1
## a[4] -1.04 0.29 -1.51 -0.61 3354 1
## a[5] -0.73 0.27 -1.17 -0.31 2449 1
## a[6] 0.23 0.27 -0.19 0.66 3304 1
## a[7] 1.82 0.39 1.20 2.44 3438 1
## bp 0.83 0.27 0.40 1.26 2193 1
## bpC -0.14 0.31 -0.63 0.33 3128 1# compare marginal posterior distributions
pairs(m10.4_map)
pairs(m10.4_stan)
10H2
# fit simpler models
m10.1 <- map(
alist(
pulled_left ~ dbinom( 1 , p ) ,
logit(p) <- a ,
a ~ dnorm(0,10)
) , data=d )
m10.2 <- map(
alist(
pulled_left ~ dbinom( 1 , p ) ,
logit(p) <- a + bp*prosoc_left ,
a ~ dnorm(0,10) ,
bp ~ dnorm(0,10)
) , data=d )
m10.3 <- map(
alist(
pulled_left ~ dbinom( 1 , p ) ,
logit(p) <- a + (bp + bpC*condition)*prosoc_left ,
a ~ dnorm(0,10) ,
bp ~ dnorm(0,10) ,
bpC ~ dnorm(0,10)
) , data=d )
# compare models with WAIC
compare(m10.1, m10.2, m10.3, m10.4_map)## WAIC pWAIC dWAIC weight SE dSE
## m10.4_map 550.5 15.8 0.0 1 18.60 NA
## m10.2 680.5 2.0 130.1 0 9.41 18.14
## m10.3 682.5 3.1 132.0 0 9.39 18.06
## m10.1 687.9 1.0 137.4 0 7.08 18.9210H3
# trace plot for stan model
plot(m10H3_stan)
# compare precis outputs
precis(m10H3_map)## Mean StdDev 5.5% 94.5%
## a 0.59 0.66 -0.47 1.65
## bp 4.24 0.90 2.81 5.67
## bv -4.59 0.96 -6.13 -3.06
## ba 1.08 0.53 0.23 1.93precis(m10H3_stan)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 0.66 0.69 -0.42 1.77 4046 1
## bp 4.64 1.01 3.06 6.18 3115 1
## bv -5.05 1.07 -6.69 -3.35 2615 1
## ba 1.14 0.54 0.27 1.98 4526 1pairs(m10H3_map)
pairs(m10H3_stan)
Overall the quadratic approximation is OK, but since the marginal posterior distributions for bp and bv are skewed, the quadratic approximation is not entirely accurate.
# the predicted probability of success and it 89% interval
mu <- link(m10H3_stan)## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]mu.mean <- apply(mu, 2, mean)
mu.HPDI <- apply(mu, 2, HPDI)
# empty plot frame
plot(0, 0, type="n", xlab="P/A/V", ylab="Probability of Success",
ylim=c(0, 1), xaxt="n", xlim=c(1,8))
axis(1, at=1:8, labels=c("1/1/1", "1/1/0", "1/0/1", "1/0/0",
"0/1/1", "0/1/0", "0/0/1", "0/0/0"))
# plot raw proportion
prop <- d$y/d$n
points(1:8, prop, pch=19)
# plot predicted probability and 89% HPDI
points(1:8-0.1, mu.mean)
points(1:8-0.1, mu.HPDI[1,], pch=3)
points(1:8-0.1, mu.HPDI[2,], pch=3)
# the predicted success count and it 89% interval
y.sim <- sim(m10H3_stan)## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]y.mean <- apply(y.sim, 2, mean)
y.HPDI <- apply(y.sim, 2, HPDI)
# empty plot frame
plot(0, 0, type="n", xlab="P/A/V", ylab="Success Count",
ylim=c(0, 30), xaxt="n", xlim=c(1,8))
axis(1, at=1:8, labels=c("1/1/1", "1/1/0", "1/0/1", "1/0/0",
"0/1/1", "0/1/0", "0/0/1", "0/0/0"))
# plot raw count
points(1:8, d$y, pch=19)
# plot predicted success count and 89% HPDI
points(1:8-0.1, y.mean)
points(1:8-0.1, y.HPDI[1,], pch=3)
points(1:8-0.1, y.HPDI[2,], pch=3)
# check trace plot
plot(m10H3_inter)
compare(m10H3_stan, m10H3_inter)## WAIC pWAIC dWAIC weight SE dSE
## m10H3_inter 93.8 4.6 0.0 0.92 12.73 NA
## m10H3_stan 98.6 4.0 4.8 0.08 13.34 4.72precis(m10H3_inter)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -0.77 1.04 -2.33 0.92 2601 1
## bp 6.58 1.43 4.42 8.87 2232 1
## bv -5.30 1.18 -7.15 -3.48 2894 1
## ba 3.44 1.24 1.51 5.42 2339 1
## bpa -2.97 1.37 -5.15 -0.81 2349 1# the predicted probability of success and it 89% interval
mu <- link(m10H3_inter)## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]mu.mean <- apply(mu, 2, mean)
mu.HPDI <- apply(mu, 2, HPDI)
# empty plot frame
plot(0, 0, type="n", xlab="P/A/V", ylab="Probability of Success",
ylim=c(0, 1), xaxt="n", xlim=c(1,8))
axis(1, at=1:8, labels=c("1/1/1", "1/1/0", "1/0/1", "1/0/0",
"0/1/1", "0/1/0", "0/0/1", "0/0/0"))
# plot raw proportion
prop <- d$y/d$n
points(1:8, prop, pch=19)
# plot predicted probability and 89% HPDI
points(1:8-0.1, mu.mean)
points(1:8-0.1, mu.HPDI[1,], pch=3)
points(1:8-0.1, mu.HPDI[2,], pch=3)
Comparing the two models with WAIC shows that the model with the interaction term gets 92% of the Akaike weights, while the model without interaction only gets 8%. So the model with the interaction term is estimated to perform much better in predicting future observations.
10H4
# check trace plot
plot(m10H4_stan)
# compare precis output
precis(m10H4_map)## Mean StdDev 5.5% 94.5%
## a -1.46 0.45 -2.18 -0.74
## bp 0.03 0.01 0.02 0.04precis(m10H4_stan)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -1.54 0.48 -2.24 -0.77 1221 1
## bp 0.03 0.01 0.02 0.04 1227 1pairs(m10H4_map)
pairs(m10H4_stan)
# plot the expected counts
PCTCOVER.seq <- seq(from=0, to=100, by=2)
pred.data <- data.frame(PCTCOVER=PCTCOVER.seq)
count.sim <- sim(m10H4_stan, data=pred.data, n=8000)## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]count.mean <- apply(count.sim, 2, mean)
count.HPDI <- apply(count.sim, 2, HPDI)
plot(SALAMAN ~ PCTCOVER, data=d, col=rangi2,
xlab="Percentage of ground cover", ylab="Counts of salamanders in each plot", pch=19)
lines(PCTCOVER.seq, count.mean)
shade(count.HPDI, PCTCOVER.seq)
# check trace plot
plot(m10H4.1)
plot(m10H4.2)
plot(m10H4.3)
compare(m10H4.1, m10H4.2, m10H4.3)## WAIC pWAIC dWAIC weight SE dSE
## m10H4.1 213.5 4.8 0.0 0.87 26.42 NA
## m10H4.2 217.6 7.6 4.1 0.11 27.27 1.44
## m10H4.3 221.9 10.4 8.4 0.01 28.18 2.42count.ensemble <-ensemble(m10H4.1, m10H4.2, m10H4.3)
count.mean <- apply(count.ensemble$sim, 2, mean)
count.HPDI <- apply(count.ensemble$sim, 2, HPDI)
observed.count <- d$SALAMAN+rnorm(n=nrow(d), mean=0, sd=0.1) # dodge points a bit
plot(count.mean ~ observed.count, xlab="Observed Count", ylab="Predicted Count",
xlim=c(0, 12), ylim=c(0, 12), col=rangi2)
abline(a=0, b=1, lty=2)
arrows(observed.count, count.HPDI[1, ], observed.count, count.HPDI[2, ],
length=0.05, angle=90, code=3)